setwd("code/")

source("insurance_source_score_1_prep.R")

#### compare expected favorability for demographics
favorlm <- lm(
    FAVOR>=3 ~ as.factor(EDUC)+scale(INCOME)+I(scale(INCOME)^2)+I(scale(INCOME)^3)+
        I(scale(AGEN))+I(scale(AGEN)^2)+I(scale(AGEN)^3)+I(scale(AGEN)^4)+I(scale(AGEN)^5)+
        BLACK+HISP+ASIAN+as.factor(ST)+MALE+RETIRED,
    data = subset(sdat, DATN < 60)
)
sdat$favor_predict <- NA
sdat$favor_predict <- predict(favorlm, newdata=sdat, type="response")


stargazer(
    favorlm,
    ## covariate.labels=c(
    ## ),
    omit=c("as.factor\\(ST\\)"),
    omit.stat=c("ser","adj.rsq", "rsq"),
    table.layout = "ts", float=F, digits=2,
    out=paste0(
        the_prefix, "/figs/",
        "tableA8_favor_demog_model.tex"
    )
)



#### plots ----
### subset data by insurance status

#### POST
sdat.1 <- subset(
    sdat, DATN > 60.5)
sdat.2 <- subset(
    sdat,
    SELFINSURE==1 & !is.na(SELFINSURE) & DATN > 60.5
)
sdat.3 <- subset(
    sdat,
    MARKET==1 & !is.na(MARKET)  & DATN > 60.5
)
sdat.4 <- subset(
    sdat,
    EMPLINSURE==1 & !is.na(EMPLINSURE) & DATN > 60.5
)
sdat.5 <- subset(
    sdat,
    MEDICAID==1 & !is.na(MEDICAID) & DATN > 60.5
)
sdat.6 <- subset(
    sdat,
    MEDICARESR==1 & !is.na(MEDICARESR) & DATN > 60.5
)
sdat.7 <- subset(
    sdat,
    OTHER==1 & !is.na(OTHER) & DATN > 60.5
)
sdat.8 <- subset(
    sdat,
    COVERED==0 & !is.na(COVERED) & DATN > 60.5
)

vec2 <- c(dim(sdat.1)[1],dim(sdat.2)[1],dim(sdat.3)[1],dim(sdat.4)[1],
          dim(sdat.5)[1],dim(sdat.6)[1],dim(sdat.7)[1],dim(sdat.8)[1])

vec2_favor <- c(
    mean(with(sdat.1, (FAVOR>=3) - favor_predict), na.rm=T),
    mean(with(sdat.2, (FAVOR>=3) - favor_predict), na.rm=T),
    mean(with(sdat.3, (FAVOR>=3) - favor_predict), na.rm=T),
    mean(with(sdat.4, (FAVOR>=3) - favor_predict), na.rm=T),
    mean(with(sdat.5, (FAVOR>=3) - favor_predict), na.rm=T),
    mean(with(sdat.6, (FAVOR>=3) - favor_predict), na.rm=T),
    mean(with(sdat.7, (FAVOR>=3) - favor_predict), na.rm=T),
    mean(with(sdat.8, (FAVOR>=3) - favor_predict), na.rm=T)
)



#### PRE
sdat.1a <- subset(
    sdat,
    DATN < 60.5
)
sdat.2a <- subset(
    sdat,
    SELFINSURE==1 & !is.na(SELFINSURE) & DATN < 60.5
)
sdat.3a <- subset(
    sdat,
    MARKET==1 & !is.na(MARKET)  & DATN < 60.5
)
sdat.4a <- subset(
    sdat,
    EMPLINSURE==1 & !is.na(EMPLINSURE) & DATN < 60.5
)
sdat.5a <- subset(
    sdat,
    MEDICAID==1 & !is.na(MEDICAID) & DATN < 60.5
)
sdat.6a <- subset(
    sdat,
    MEDICARESR==1 & !is.na(MEDICARESR) & DATN < 60.5
)
## sdat.6a <- subset(
##     sdat,
##     AGEN>=65 & !is.na(AGEN) & !is.na(MEDICARESR) & DATN < 60.5
## )
sdat.7a <- subset(
    sdat,
    OTHER==1 & !is.na(OTHER) & DATN < 60.5
)
sdat.8a <- subset(
    sdat,
    COVERED==0 & !is.na(COVERED) & DATN < 60.5
)

vec1 <- c(dim(sdat.1a)[1],dim(sdat.2a)[1],dim(sdat.3a)[1],dim(sdat.4a)[1],
          dim(sdat.5a)[1],dim(sdat.6a)[1],dim(sdat.7a)[1],dim(sdat.8a)[1])

vec1_favor <- c(
    mean(with(sdat.1a, (FAVOR>=3) - favor_predict), na.rm=T),
    mean(with(sdat.2a, (FAVOR>=3) - favor_predict), na.rm=T),
    mean(with(sdat.3a, (FAVOR>=3) - favor_predict), na.rm=T),
    mean(with(sdat.4a, (FAVOR>=3) - favor_predict), na.rm=T),
    mean(with(sdat.5a, (FAVOR>=3) - favor_predict), na.rm=T),
    mean(with(sdat.6a, (FAVOR>=3) - favor_predict), na.rm=T),
    mean(with(sdat.7a, (FAVOR>=3) - favor_predict), na.rm=T),
    mean(with(sdat.8a, (FAVOR>=3) - favor_predict), na.rm=T)
)

mat <- cbind(vec1[2:8]/vec1[1],vec2[2:8]/vec2[1])
rownames(mat) <- c(
    "Self-Purchased","Exchange","Employer-Insured",
    "Medicaid","Medicare","Other","Uninsured"
)

mat_favor <- cbind(vec1_favor[2:8],vec2_favor[2:8])
rownames(mat_favor) <- c(
    "Self-Purchased","Exchange","Employer-Insured",
    "Medicaid","Medicare","Other","Uninsured"
)

## mat2 <- mat[order(mat[,1],decreasing=T),]
mat2 <- mat[order(mat[,1],decreasing=F),]
bp <- barplot(t(mat2),beside=T)

pdf(paste0(the_prefix, "/figs/figureA1_kff-insurance-by-type-20190114.pdf"),width=7, height=4.5)
par(mar=c(10,3,2,2))
x <- barplot(
    t(mat2),
    beside=T, cex.axis=1.5, las=2, cex.names=1.25,
    xaxt="n", col=c(gray(0.9),gray(0.3))
)
text(x=x[1,]-.25, y=-0.15, rownames(mat2), xpd=TRUE, srt=45, cex=1.25)
##
text(y=.25,x=1.5,cex=1.25,"Before \n1/14")
text(y=.25,x=4,cex=1.25,"After \n1/14")
##
arrows(y0=.2,x0=1.5,x1=bp[1,1],y1=.05)
arrows(y0=.2,x0=4,x1=bp[2,1],y1=.05)
dev.off()


pdf(paste0(the_prefix, "/figs/figuresA3_A4_kff-insurance-favor-by-type-20190114.pdf"),width=7, height=4.5)
## favorability
mat_favor2 <- mat_favor[order(mat[,1],decreasing=F),]
##
par(mar=c(8,8,2,2))
x_favor <- barplot(
    t(mat_favor2),
    beside=T, cex.axis=1.5, las=2, cex.names=1.25,
    xaxt="n", col=c(gray(0.9),gray(0.3))
)
text(x=x[1,]-.25, y=-0.18, rownames(mat_favor2), xpd=TRUE, srt=45, cex=1.25)
mtext("Actual - Expected\nACA Favorability", side=2, padj=-1.75, cex=1.5)
legend(
    "topright", bty="n", legend=c("Before 1/14", "After 1/14"),
    col=c(gray(0.9),gray(0.3)), pch=15, pt.cex=2, cex=1.5
)
abline(h=0, col="red", lwd=2)
##
## favorability x size
mat_favor2_alt <- (mat2 * mat_favor2)
par(mar=c(8,9,2,2))
x_favor <- barplot(
    t(mat_favor2_alt),
    beside=T, cex.axis=1.5, las=2, cex.names=1.25,
    xaxt="n", col=c(gray(0.9),gray(0.3))
)
text(
    x=x[1,]-.25, y=-0.015, rownames(mat_favor2_alt),
    xpd=TRUE, srt=45, cex=1.25
)
mtext(
    "Actual - Expected ACA\nFavorability x Pop Size",
    side=2, padj=-2.25, cex=1.5
)
legend(
    "bottomleft", bty="n", legend=c("Before 1/14", "After 1/14"),
    col=c(gray(0.9),gray(0.3)), pch=15, pt.cex=2, cex=1.5
)
abline(h=0, col="red", lwd=2)
dev.off()


#### Table A1
n <- 8
rmat <- matrix(NA,n,8)
for(i in 1:n){
	txt <- paste("dta.sub <- sdat.",i,sep="")
	eval(parse(text=txt))
	rmat[i,1] <- mean(dta.sub$INCOME,na.rm=T)
	rmat[i,2] <- mean(dta.sub$EDUC,na.rm=T)
	rmat[i,3] <- mean(dta.sub$AGEN,na.rm=T)
	rmat[i,4] <- mean(dta.sub$MALE,na.rm=T)
	## rmat[i,5] <- mean(dta.sub$ASIAN,na.rm=T)
	rmat[i,5] <- mean(dta.sub$BLACK,na.rm=T)
	rmat[i,6] <- mean(dta.sub$HISP,na.rm=T)
	rmat[i,7] <- mean(dta.sub$PID5,na.rm=T)
	rmat[i,8] <- dim(dta.sub)[1]
}
rownames(rmat) <- c("All","Own Ins.","Used Exchanges","Employer-Based Ins.",
	"Medicaid","Medicare","Other","Uninsured")
colnames(rmat) <- c(
    "Income","Educ.","Age","Male","Black","Hispanic","Party ID","N"
)

## round(rmat,digits=3)
xtable(rmat,digits=c(0,1,1,0,rep(2,3),1,0))
print(
    xtable(rmat,digits=c(0,1,1,0,rep(2,3),1,0)),
           file=paste0(the_prefix, "/figs/tableA1_covariates_post_implementation.tex")
    )
